function [h1,l2, h2,L2,H1]=bberrorb(dofs,mesh,x,desc)
% if nargin<=2
%    ux=[];
%    uy=[];
% end

%comput the total striff matrix

[dofs n_dof]=sort_dof_dis2(mesh.T,size(mesh.T,1),mesh.d);
dcurrent=mesh.d;
dmax=max(dcurrent);
    for din=1:dmax
Mat{din} = build(dmax-1);
[cr,desc]=get_auxillary_mat(din);
Id = diag(ones((din+1)*(din+2)/2,1));
Du_cell{din} = din*de_cast_step(Id,din,1,0,-1,desc); % the direction derivate on v3-v1
Dv_cell{din} = din*de_cast_step(Id,din,0,1,-1,desc);  % the direction derivate on v3-v2
Duu_cell{din}= (din-1)*de_cast_step(Du_cell{din},din-1,1,0,-1,desc);
Duv_cell{din}= (din-1)*de_cast_step(Du_cell{din},din-1,0,1,-1,desc);
Dvv_cell{din}= (din-1)*de_cast_step(Dv_cell{din},din-1,0,1,-1,desc);
    end
[qp,qw,b1,b2,b3]=quadrule;

% reserve spaces for entry
T=mesh.T;
V=mesh.V;
posT=size(T,1);
pos = 1;
d_plot=3;
m_plot = (d_plot+1)*(d_plot+2)/2;
for k = 1:posT
    d=dcurrent(k);
    Du=Du_cell{d};
    Dv=Dv_cell{d};
     Duu=Duu_cell{d};
     Duv=Duv_cell{d};
     Dvv=Dvv_cell{d};
    
    V1=V(T(k,1),:);V2=V(T(k,2),:);V3=V(T(k,3),:);
    x13 = V1(1)-V3(1);
    y13 = V1(2)-V3(2);
    x23 = V2(1)-V3(1);
    y23 = V2(2)-V3(2);
    J = x13*y23 - x23*y13; 
     qp_glob = zeros(size(qp));
%       qp_glob(:,1) = V3(1) + x13*qp(:,1) + x23*qp(:,2);
%       qp_glob(:,2) = V3(2) + y13*qp(:,1) + y23*qp(:,2);  
      % evaluate at integrate points, 
        qp_glob(:,1) = b1*V1(1) + b2*V2(1) + b3*V3(1);
      qp_glob(:,2) =  b1*V1(2) + b2*V2(2) + b3*V3(2); 
%       fvec = feval(u, qp_glob(:,1), qp_glob(:,2),4);
%        fxvec = feval(ux, qp_glob(:,1), qp_glob(:,2));
%         fyvec = feval(uy, qp_glob(:,1), qp_glob(:,2));
%          fxxvec = feval(uxx, qp_glob(:,1), qp_glob(:,2));
%         fxyvec = feval(uxy, qp_glob(:,1), qp_glob(:,2));
%         fyyvec=feval(uyy, qp_glob(:,1), qp_glob(:,2));
    loc_dof = get_tri_dof(dofs,k);
   loc_c=mesh.c(loc_dof);
    fbb=vdm23(d,b1,b2,b3)*loc_c;
     fbbx=(vdm23(d-1,b1,b2,b3)*Du)*loc_c*y23/J-(vdm23(d-1,b1,b2,b3)*Dv)*loc_c*y13/J;
     fbby=-(vdm23(d-1,b1,b2,b3)*Du)*loc_c*x23/J+(vdm23(d-1,b1,b2,b3)*Dv)*loc_c*x13/J;
    fbbxx=vdm23(d-2,b1,b2,b3)*Duu*loc_c*(y23/J)^2-2*vdm23(d-2,b1,b2,b3)*Duv*loc_c*(y23*y13)/J^2 +vdm23(d-2,b1,b2,b3)*Dvv*loc_c*(y13/J)^2 ;
     fbbyy=vdm23(d-2,b1,b2,b3)*Duu*loc_c*(x23/J)^2-2*vdm23(d-2,b1,b2,b3)*Duv*loc_c*x23*x13/J/J +vdm23(d-2,b1,b2,b3)*Dvv*loc_c*(x13/J)^2;
     fbbxy=vdm23(d-2,b1,b2,b3)*Duu*loc_c*(-y23*x23)/J/J+vdm23(d-2,b1,b2,b3)*Duv*loc_c*(x23*y13+y23*x13)/J^2 +vdm23(d-2,b1,b2,b3)*Dvv*loc_c*(-x13*y13)/J^2;
     fbbxx=vdm23(d-2,b1,b2,b3)*Duu*loc_c*(y23/J)^2-2*vdm23(d-2,b1,b2,b3)*Duv*loc_c*(y23*y13)/J^2 +vdm23(d-2,b1,b2,b3)*Dvv*loc_c*(y13/J)^2 ;
   loc_c=x(loc_dof);
     fvec=vdm23(d,b1,b2,b3)*loc_c;
     fxvec=(vdm23(d-1,b1,b2,b3)*Du)*loc_c*y23/J-(vdm23(d-1,b1,b2,b3)*Dv)*loc_c*y13/J;
     fyvec=-(vdm23(d-1,b1,b2,b3)*Du)*loc_c*x23/J+(vdm23(d-1,b1,b2,b3)*Dv)*loc_c*x13/J;
    fxxvec=vdm23(d-2,b1,b2,b3)*Duu*loc_c*(y23/J)^2-2*vdm23(d-2,b1,b2,b3)*Duv*loc_c*(y23*y13)/J^2 +vdm23(d-2,b1,b2,b3)*Dvv*loc_c*(y13/J)^2 ;
     fyyvec=vdm23(d-2,b1,b2,b3)*Duu*loc_c*(x23/J)^2-2*vdm23(d-2,b1,b2,b3)*Duv*loc_c*x23*x13/J/J +vdm23(d-2,b1,b2,b3)*Dvv*loc_c*(x13/J)^2;
     fxyvec=vdm23(d-2,b1,b2,b3)*Duu*loc_c*(-y23*x23)/J/J+vdm23(d-2,b1,b2,b3)*Duv*loc_c*(x23*y13+y23*x13)/J^2 +vdm23(d-2,b1,b2,b3)*Dvv*loc_c*(-x13*y13)/J^2;
     fxxvec=vdm23(d-2,b1,b2,b3)*Duu*loc_c*(y23/J)^2-2*vdm23(d-2,b1,b2,b3)*Duv*loc_c*(y23*y13)/J^2 +vdm23(d-2,b1,b2,b3)*Dvv*loc_c*(y13/J)^2 ;
     
     
     L2(k)=((fbb-fvec).^2)'*qw*J;
    H1(k)= ((fbby-fyvec).^2)'*qw*J+((fbbx-fxvec).^2)'*qw*J;
    H2(k)=((fbbxx-fxxvec).^2)'*qw*J+((fbbxy-fxyvec).^2)'*qw*J+((fbbyy-fyyvec).^2)'*qw*J;%
end
h1=sqrt(sum(L2)+sum(H1));
l2=sqrt(sum(L2));
h2=sqrt(sum(L2)+sum(H1)+sum(H2));
L2;
H1;
